Skip to content

Round to_f64/to_f32 once, to the subnormal-aware precision - #91

Merged
cmpute merged 3 commits into
cmpute:masterfrom
gaoflow:fix-subnormal-double-rounding
Jul 22, 2026
Merged

Round to_f64/to_f32 once, to the subnormal-aware precision#91
cmpute merged 3 commits into
cmpute:masterfrom
gaoflow:fix-subnormal-double-rounding

Conversation

@gaoflow

@gaoflow gaoflow commented Jul 18, 2026

Copy link
Copy Markdown
Contributor

FBig/DBig → float conversion rounds the source to a fixed 53/24-bit intermediate and then lets f{32,64}::encode round again into the format's exponent-dependent subnormal grid. That second rounding is wrong in two ways. This is the subnormal counterpart of the normal-range double-rounding that was recently addressed.

Subnormal double-rounding (1 ULP). A subnormal has fewer than 53/24 significant bits — its spacing is fixed at 2^-1074 / 2^-149 — so the fixed intermediate keeps bits below the subnormal grid. For a value just past a subnormal halfway, the intermediate rounds to the halfway (dropping the excess that lifts it above) and encode then rounds half-to-even down, one ULP below the correctly-rounded value:

// value = (2m+1)·2^-1075 + 2^-1135, just above the halfway between subnormals m and m+1
let sig = (IBig::from(2 * m + 1) << 60) + IBig::ONE;
let v = FBig::<HalfEven, 2>::from_parts(sig, -1135);
v.to_f64(); // returned m·2^-1074 (rounded to even); correct is (m+1)·2^-1074

Debug panic / release double-round on high-precision inputs. to_f64/to_f32 feed the unrounded source significand into the base-changing division, whose debug_assert!(lhs.digits() <= self.precision + rhs.digits()) fires whenever the significand is wider than the target precision:

DBig::from_str("123456789012345678.9012345678901")?.to_f64(); // panics in debug builds

In release the assert is compiled out and the oversized quotient is re-rounded downstream — a silent double rounding.

Fix. Round the source once, straight to the target's precision at its own magnitude (min(53, e + 1075) bits for f64 at binary exponent e, analogously for f32), so encode re-rounds nothing. To keep that single rounding correct without an exact arbitrary-precision conversion, the source is first taken to a fixed, generous width with round-to-odd; rounding that down to the final width then reproduces the correctly-rounded result for every rounding mode. The division now rounds an over-wide quotient down to the requested precision instead of assuming the caller pre-bounded the dividend, so the assertion is removed and high-precision inputs convert without panicking (and without the release double-round).

Correctness holds to within the base conversion's working precision; a value closer to a boundary than the conversion resolves stays a pre-existing limitation of finite-precision base changing (shared with with_base). The base-2 / Repr paths are exact for all inputs.

Tests cover subnormal halfways nudged just past the boundary across the exponent range, in both base 2 and base 10, for f64 (oracle: float()) and f32 (oracle: exact Fraction → nearest 2^-149, not numpy.float32, which itself double-rounds), plus the high-precision decimals that previously panicked. cargo test --workspace --exclude dashu-python, cargo clippy --all-features --all-targets --workspace --exclude dashu-python -- -D warnings, and cargo fmt --all -- --check all pass.

The conversions rounded the source to a fixed 53/24-bit intermediate, then
the bit encoding rounded again into the exponent-dependent subnormal grid.
A value just past a subnormal halfway lands on the halfway after the first
rounding and the second rounds to even -- one ULP low. Round the source
once, straight to the target's precision at its magnitude (fewer than 53/24
bits for subnormals) through a round-to-odd base conversion, so the single
rounding is mode-correct and never double-rounds.

Also round an over-wide quotient down to the requested precision in the
base-changing division instead of requiring the caller to pre-bound the
dividend: high-precision decimals (e.g. "123456789012345678.9012345678901")
previously panicked in debug and double-rounded in release when converted.
- Cover negative inputs in the to_f{32,64} subnormal-halfway and
  high-precision tests (sign bit flip); the new over-wide-quotient path
  and subnormal rounding are sign-aware but were only tested positive.
- Reword the convert_base_odd doc: rounding down to width-2 reproduces
  the correctly-rounded value for every rounding mode (not just nearest).
- Fold f64/f32_significand_bits into a shared significand_bits helper
  parameterized by max_bits and the subnormal exponent.

Co-Authored-By: Claude <noreply@anthropic.com>
@cmpute
cmpute force-pushed the fix-subnormal-double-rounding branch from 8416569 to c8b79e1 Compare July 22, 2026 15:14
repr_div carries a documented precondition — lhs.digits() <= precision +
rhs.digits() — that callers must uphold (Context::div already pre-shrinks
the dividend to exactly this bound). The previous fix instead relaxed the
contract: it dropped the debug_assert and taught repr_div itself to round
an over-wide quotient down. That worked but departed from the layering —
repr_div is a lightweight hot kernel whose complexity belongs in the
informed caller, and the assert loss demoted a real invariant to an
implicit convention.

Restore repr_div's contract and assert, and fix the actual offender: the
to_f64/to_f32 path reaches repr_div via convert_base's small-exponent
division, which fed an oversized significand straight in. convert_base now
pre-shrinks the dividend to den.digits() + precision before dividing,
mirroring Context::div. Behavior is unchanged (the high-precision and
subnormal tests still pass, now through the caller path, with the
debug_assert active).

Co-Authored-By: Claude <noreply@anthropic.com>
@cmpute
cmpute merged commit 4b0bd88 into cmpute:master Jul 22, 2026
@cmpute

cmpute commented Jul 22, 2026

Copy link
Copy Markdown
Owner

Thanks! Merged with some personal modifications

cmpute added a commit that referenced this pull request Jul 24, 2026
* Address opendp-num fuzz findings (DASHU-007/008/015/022/026)

Reviewed the opendp-num differential/property-fuzz findings against the
current source. DASHU-020/021 are already fixed by PR #91; DASHU-023/024
(directed rounding at FBig's exponent-range extremes) are tracked with
TODO comments and deferred. This commit fixes the remaining five.

DASHU-015 (dashu-int, incorrect-result): `to_f64_small` reported an inexact
conversion as `Approximation::Exact` at the `DoubleWord::MAX` boundary. The
exactness test used a saturating `f as DoubleWord` round-trip, which clamps
back to `DoubleWord::MAX` when the value rounds up to `2^BITS`. Detect that
saturation case before the comparison; apply the same guard to the 16/32-bit
RefLarge -> f64 fast path.

DASHU-026 (dashu-float, debug panic/OOM): `round_fract`'s debug assert built
`B^precision`; for a sparse sticky tail (precision == the exponent gap, e.g.
from `exp_m1` of a large-magnitude input) this exhausted memory in debug
builds. Replace it with a `log2_bounds`-based precondition check that
allocates nothing and only fires on a proven violation.

DASHU-022 (dashu-float, panic): `assert_limited_precision` fired before the
exact zero/one shortcuts, so values carrying unlimited precision (precision
0) -- `FBig::try_from(0.0)`, `FBig::ONE`/`ZERO` -- panicked in exp/exp_m1/
sqrt/ln/ln_1p despite mathematically exact results. Hoist the exact
shortcuts above the assertion.

DASHU-008 (dashu-int, panic): GCD of two large, similarly-sized integers
aborted with "internal error: not enough memory allocated". The scratchpad
was reserved once from the initial operand lengths, but each euclidean
step's division dispatches on the current lengths, so a later lopsided
Burnikel-Ziegler step (divisor > 48 words, quotient > 32 words) was
under-reserved. Size it for the worst case reachable in the loop:
`mul::memory_requirement_up_to(rhs_len, rhs_len / 2)`. Applies to both the
gcd and extended-gcd reservations.

DASHU-007 (dashu-float, missing feature): add a correctly-rounded
`FBig::log2` / `Context::log2` / `CachedFBig::log2`, computed as
`ln(x)/ln(2)` at an elevated working precision. Previously only the
f32-precision `log2_bounds` estimate existed, so directed log2 was wrong by
many ULPs. The result magnitude tracks the division's error amplification,
so a few guard digits certify the final round across the whole range.

Each fix includes a regression test; the DASHU-008 test reproduces the exact
memory.rs panic under the old code.

Co-Authored-By: Claude <noreply@anthropic.com>

* Drop opaque finding IDs from code/test comments

The DASHU-NNN tokens reference an external fuzz-findings corpus that other
developers don't have context for. Rephrase the affected comments (two
production TODOs and several test docstrings) to describe the limitation or
regression in self-contained terms.

Co-Authored-By: Claude <noreply@anthropic.com>

* Use fixed operands in GCD tests; ban random generators in unit tests

Replace the xorshift-based `big_from_seed` in the two DASHU-008 regression
tests (`test_gcd_large_lopsided_reduction`, `rbig_reduce_large_lopsided`)
with fixed all-ones operands sized by word count. The tests still reproduce
the exact `memory.rs:150` panic under the old code (verified) but now use
deterministic inputs.

Add a Code-style rule to AGENTS.md: in-crate tests must use fixed,
deterministic inputs — never random/property-test generators — except when
testing the `rand` integration itself. Randomized input belongs under
`fuzz/`.

Co-Authored-By: Claude <noreply@anthropic.com>

* Gate exp_m1 OOM regression test to 64-bit

On 32-bit targets `isize` tops out at ~2.1e9, so for x = -2^62 the value
floor(x/ln2) ≈ -6.6e18 overflows isize and exp_internal takes its (deferred)
overflow branch — returning Exact(-1) — instead of the round_fract path this
test targets. A sharp OOM needs an exponent gap large enough that 2^gap
exceeds memory, yet still fitting isize; that window only exists on 64-bit.
Gate the test accordingly. The underlying round_fract fix is arch-independent.

Co-Authored-By: Claude <noreply@anthropic.com>

* Restore TODO markers on the deferred powi/exp limitation comments

The earlier ID cleanup dropped the TODO tag along with the issue ID.
Keep `TODO:` so the limitations still surface in a TODO grep.

Co-Authored-By: Claude <noreply@anthropic.com>

---------

Co-authored-by: Jacob Zhong <jacob@rimbot.com>
Co-authored-by: Claude <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants